Nonlinearity and Multifractality of Climate Change in the Past 420,000 Years 
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Evidence of past climate variations are stored in ice and indicate glacial-intergiacial cycles charac- 
terized by three dominant time periods of 20kyr, 40kyr, and lOOkyr. We study the scaling properties 
of temperature proxy records of four ice cores from Antarctica and Greenland. These series are long- 
range correlated in the time scales of l-100kyr. We show that these series are nonlinear as expressed 
by volatility correlations and a broad multifractal spectrum. We present a stochastic model that 
captures the scaling and the nonlinear properties observed in the data. 



PACS numbers: 92.70.Gt, 05.40.-a, 92.40. Cy 

Abundant geological evidence indicates that temper- 
atures varied from the cold of ice ages to the warmth 
of interglacial periods. In the last 800,000 years 
(800kyr) there is strong evidence for a dominant glacial- 
intergiacial cycle of lOOkyr, with weaker secondary cy- 
cles of 40kyr and 20kyr Each lOOkyr cycle consists of 
gradual cooling for ~ 90kyr followed by rapid warming 
during ~ lOkyr. "Milankovitch forcing" , which refers to 
changes in insolation (solar radiation) due to variations 
in the precession, obliquity (tilting), and eccentricity of 
Earth's orbit g] are thought to play an important role in 
glacial dynamics. These orbital variations are character- 
ized by periods of 20kyr, 40kyr, and lOOkyr, respectively. 
The 20kyr and 40kyr periods in the climate records are 
generally believed to be a linear response of the climate 
system to insolation variations. In contrast, the weak- 
ness of the variations in solar radiation at the lOOkyr 
timescale has lead to the generally accepted conclusion 
that the glacial-intergiacial oscillations at this timescale 
are most likely not a direct linear response of the climate 
system to external solar variations || . 

Many deterministic theories have been developed to 
explain the glacial-intergiacial lOOkyr variability; the ma- 
jority suggest that the lOOkyr period is a result of self- 
sustained nonlinear mechanisms (see, e.g., [§-[|]). Other 
studies proposed that climate variations are stochas- 
tic and follow scaling laws — the Milankovitch periods 
are second-order perturbations (e.g. ||-f|). Importantly, 
both the deterministic and stochastic mechanisms still 
assume that the variations on time scales below lOOkyr 
down to lOkyr are linear. The objectives of the present 
study are to quantify the degree of nonlinearity of climate 
dynamics within the time scales of l-100kyr and to pro- 
vide statistical characteristics of the proxy records which 
can serve as a test for distinguishing between existing 
climate models ||. 

We study the correlation (scaling) properties of climate 
records of the past 420kyr. We show that temperature 
variations are long-range correlated suggesting that the 



Milankovitch periods are indeed secondary and (contrary 
to common belief [||) that climate dynamics of all time 
scales below lOOkyr down to lkyr are highly nonlinear. 
In addition, we quantify the degree of nonlinearity in the 
climate records and suggest a possible stochastic nonlin- 
ear mechanism for our findings. 

Our analysis is based on isotope records obtained from 
four ice cores, Vostok and Taylor Dome from Antarc- 
tica, and GISP (Greenland-Ice-Sheet-Project) and GRIP 
(Greenland-Ice-Project) from Greenland ||. Measure- 
ments of oxygen and hydrogen isotope ratios (S 18 and 
5D) of the ice at different depths in the core provide a 
proxy record of temperature p0[ when the ice was formed 
(Fig. |l|a). These records extend back to 100-420kyr. 

Fourier analysis is the standard method for studying 
long-range correlations in time series. When the power 
spectrum follows scaling laws, S(f) ~ l//' 3 Q3 > 0), the 
series is long-range correlated [ p"l[ . However, the power 
spectrum might yield an inaccurate estimation of the 
scaling exponent due to constant or polynomial trends 
that are not necessarily related to the intrinsic dynamics 
|fl2f . We therefore use the detrended fluctuation analysis 
(DFA) Jl2| ; the mth order DFA eliminates polynomial 
trends of order m — 1 from the data and provides a more 
accurate estimation of the scaling exponent Jl^,|l^]. If 
the root mean square fluctuation function, F(n), is pro- 
portional to n a , where n is the window scale, the series 
is long range correlated (/3 = 2a — 1). For a random se- 
ries a = 0.5 while for correlated (or anticorrelated) series 
a > 0.5 (or a < 0.5). We begin our analysis with the Vos- 
tok ice core and find that temperature changes are highly 
correlated in the time range l-100kyr with a scaling ex- 
ponent a sa 1.5 (Fig. ^a), consistent with the previously 
reported power spectrum exponent (3 — 2 ||-||] . 

Next, we analyze the nonlinear properties of the ice 
core record. We define a process to be linear if it is pos- 
sible to reproduce its statistical properties (such as the 
third moment) from the power spectrum and the proba- 
bility distribution alone, regardless of the Fourier phases 
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p5[ . This definition includes autoregression processes 

(x n = J2iti a i x n-i + Y%=o ^iVn-i where i] is Gaussian 
white noise) and fractional Brownian motion; the output, 
x n , of these processes may undergo monotonic nonlinear 
transformations s n = s(x n ) and still be linear. Processes 
which are not linear are defined as nonlinear |fl6|| . 

Long-range correlations in the temperature time se- 
ries, Tj, reflect linear aspects of Tj. Long-range cor- 
relations in the magnitudes of temperature increments, 
|ATj| = \T l+1 -Ti\ (Fig. |b), which we define as volatility, 
indicate nonlinearity of the underlying process [ ]l6| , |l7| |. 
Linear series have uncorrelated |ATi| series while nonlin- 
ear processes that follow a scaling law exhibit long-range 
correlations in the magnitude series |ATj|. We find that 
the magnitude series | AX!;| is highly long-range correlated 
within the time range l-100kyr (Fig. |b). Thus, the un- 
derlying process is nonlinear |l7|] . The value of the cor- 
relation exponent quantifies the degree of nonlinearity in 
the ice core record. Correlations in the magnitude series 
indicates that the magnitude series is "clustered", i.e., 
large magnitude is more likely to be followed by a large 
magnitude, as can be seen in Fig. |l|b. These clusters 
may be associated with abrupt warming events known as 
Dansgaard-Oeschger events Q . 

To demonstrate that the correlations in the magnitude 
series are related to the nonlinearity of the underlying 
process we apply a surrogate data test for nonlinearity 
that preserves both the power spectrum and the his- 
togram of the temperature increment series AXi |ll| . The 
surrogate series has random Fourier phases; the nonlin- 
earities that are stored in the phases are destroyed. We 
find that the magnitude series obtained from the surro- 
gate series is indeed uncorrelated (Fig. ^o) confirming 
that the original series is nonlinear within l-100kyr. 

Correlations in the magnitude series |ATi| can be re- 
lated to the width of the multifractal spectrum p||l7j ]. 
We calculate the exponents r(q) of different moments q 
for the ice core data and find that r(q) is a nonlinear func- 
tion of q (Fig. ||c) , indicating that the temperature series 
is multifractal. We also perform multifractal analysis on 
the surrogate data and find that its r(q) is almost linear. 
The multifractal spectrum, D(h) (Fig. ||d), is broad for 
the original data and narrower for the surrogate data. 
The broadness of the multifractal spectrum may also be 
used to quantify the degree of multifractality, and thus 
the degree of nonlinearity, in the data fl7|| . 

We repeat the above analysis for the other three ice 
cores and obtain similar results (Table ||). Although the 
DFA exponents of the original series are smaller for the 
Greenland cores, the magnitude series exponents are al- 
most the same for all cores. We thus conclude that cli- 
mate dynamics is nonlinear for time scales of few thou- 
sands of years up to lOOkyr. 

To understand the mechanism that may contribute 
to the nonlinearities observed in the data we modified 



a model for ice-volume evolution recently suggested by 
Wunsch Ice volume V was observed to be nega- 
tively correlated with temperature T ~ — y |jj and thus 
a model for ice accumulation may serve, indirectly, as 
a model for temperature dynamics. The Wunsch model 
can be summarized as follows: The ice-sheet builds ran- 
domly up to a critical volume where it breaks up rapidly. 
Then, growth begins again. By construction, this model 
is a random walk up to a time scale corresponding to the 
critical ice volume, followed by a crossover to random be- 
havior for larger time scales. However, this model does 
not reproduce the nonlinearity in the data as defined and 
found above. 

The assumptions of our model (Fig. |§|) are: 

(i) The ice volume V changes with steps 6 + b/V; i.e., 
V(t + dt) = V{t)+5{t)+b/V(t). 

(ii) When ice volume V "crosses" a critical volume £, b 
is set to be negative, b = 62 < 0. Ice volume is considered 
to lie between 0.01£ and £. When V = 0.01£, b becomes 
b = b\ > till V exceeds the threshold £. 

(iii) The ice accumulation increments S are the product 
of two stochastic inputs, 5(t) = CifyVitt)- C an d V are 
Gaussian distributed random variables with zero mean 
and unit variance. 

(iv) Random switching between the states rj^s is con- 
trolled by i(t) which is equal to [l(t)] where [•] stands for 
the closest integer value. I is a random walk described 
by, l(t + dt) = l(t) + Cui(t), where C is the switching 
range and lu is another Gaussian random variable with 
zero mean and unit variance (see Fig. |^a). 

Assumptions (i) and (ii) describe the random growth 
(with b — bi > 0) of the ice-sheet and its rapid breakup 
(with b = 62 < 0) after crossing the critical volume £. 
The 1/V term in assumption (i) mimics the reduced ice 
accumulation for large ice- volume |Q . Ice volume changes 
5 result from two interacting random inputs [assumption 
(iii)] where one, £, may represent the atmosphere, the 
other, 77, the ocean, and the product, £77, the atmosphere- 
ocean interaction. 

In our simulations of the model we use the following 
values, dt = O.lkyr, £ = 1.5y/90kyr/dt, bi — 1, 62 = —3, 
and C = 0.27. We choose the value 90kyr in £ so that 
on average V will grow from zero to £ after 90kyr/dt 
steps. We choose the values of 61 and &2 such that the 
ice-sheet grows slowly and breakup rapidly. The values of 
£, b\, and 62 are constrained by the natural record. The 
switching parameter C determines the number of states 
X] for a given number of steps; the number of different 
states is proportional to the square root of the number 
of steps (e.g., ~ 3 states for 4kyr and ~ 5 for lOkyr). 

An example of an arbitrary 400kyr time series obtained 
by the model is shown in Fig. ||b. The scaling of the 
model's V series (Fig. |]a) indicates random walk behav- 
ior with exponent a — 1.5 (as for the Vostok core). The 
magnitude series |AV| is highly correlated with exponent 
- 0.8 (Fig. |b) as for the ice core data (Fig. |^b and Ta- 
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ble |) . The surrogate data test applied to the A V series 
changes the magnitude series into an uncorrelated one, 
indicating the nonlinearity of the model. This nonlin- 
earity is mainly due to the product of the inputs rj, £ 
in assumption (iii). The multifractal spectrum is broad 
where, as with the ice core data (Fig. ||c,d), the expo- 
nents for negative moments, r(q < 0), mainly contribute 
to its broadness (Fig. f|c,d). After the surrogate data 
test the series becomes linear and statistically different 
from the original data. 

This simple model reproduces the statistical character- 
istics of the ice core data under consideration. Although 
the natural system is undoubtedly more complex, we con- 
jecture that the model variables may be associated with 
specific aspects of Earth's climate system although our 
model cannot uniquely identify them. One of the ran- 
dom inputs, 77, thus may represent the influence of the 
deep ocean on ice accumulation since the state of the 
deep ocean is known to have impact on glaciation (e.g., 
^|). The other random input, C, may represent the net 
atmospheric influence affecting ice accumulation (result- 
ing from, e.g., variations in eddy transport, cloudiness, 
incoming solar radiation, ablation). We assume that the 
ocean has several states with a tendency to return to 
previous states, as does rj in the model; a possible ex- 
ample for such "switching" mechanism is the deep ocean 
circulation which has few states with possible switching 
between them (e.g. Q). 

We conclude that climate changes in the time range of 
l-100kyr are long-range correlated confirming the major 
role of stochasticity in climate Moreover, our re- 

sults suggest that the underlying dynamics in the time 
scales of l-100kyr is nonlinear. This nonlinearity is spec- 
ified and quantified by strong long-range correlations in 
the magnitudes of temperature changes and in a broad 
multifractal spectrum. Our simple stochastic model sug- 
gests that the nonlinearity can be the result of only two 
random processes that interact with each other. Glacia- 
tion models may be generally categorized into two main 
alternatives: (i) linear mechanisms that are driven by 
stochastic forcing (e.g. |@,||), and (ii) nonlinear mecha- 
nisms without stochastic forcing (e.g. ||[|]). Our results 
and model suggests a third alternative — nonlinear mech- 
anism that inherently involves stochastic forcing. Our re- 
sults raise a new challenge for the many climate models, 
and may help guide development of better climate mod- 
els, which include both periodic and stochastic elements 
of climate change. 
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FIG. 1. (a) Isotopic temperature record from the Vostok 
ice core jj]; the temperature T is calculated from the hydrogen 
isotope ratio, (b) A typical example of the magnitude of 
temperature changes | ATi | . The magnitude series is clustered 
— big magnitudes are likely to be followed by big magnitudes 
suggesting the presence of correlations in the | ATi | series. 
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FIG. 2. (a) The root mean square fluctuation -F(n) as a 
function of window scale n in kyr for the Vostok temperature 
proxy data indicates strong correlations (■). Before applying 
the second order DFA we sampled the data at O.lkyr inter- 
vals |f| . The surrogate series (gray triangles) exhibits almost 
identical scaling confirming that correlations in the T series 
are a linear measure, (b) F(n) for the magnitude series, \AT\, 
indicates strong correlations (■) . The magnitude series of 
the surrogate data (gray triangles) is uncorrelated (with ex- 
ponent 0.5) demonstrating the nonlinearity of the data, (c) 
The multifractal analysis uses the wavelet transform modu- 
lus maxima metho d hs| , with the 8-tap Daubechies discrete 
wavelet transform [EOjTThe exponents r(q) are estimated by 



scaling function Z q (n) 



Mi) 



m (0.8kyr < n <25.6kyr). 



The curvature in r(q) reflects the multifractality of the tem- 
perature series (•). The r(q) of the surrogate series (o) is 
linear (10 realizations; the average ± 1 standard deviation is 
shown), (d) The multifractal spectrum, D(h) = hq — r(q) 
(h = dr/dq), is much broader for the original data (•) com- 
pared to the average D(h) of the surrogate data (o) confirming 
that the underlying dynamics is nonlinear. 



TABLE I. Scaling 
(O.lkyr sampling ll 31 



results of the cores under consideration 
). The DFA exponents for the original 



(or) and magnitude (c*| at| ) series are obtained for window 
scales between lkyr and ~ 1/4 of the series length. The mul- 
tifractal analysis is not accurate for short series, GISP, GRIP, 
Taylor Dome, and is not presented; nonetheless, the multifrac- 
tal spectrum for these cores is broad as for Vostok (Fig. 0). 
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FIG. 3. (a) An illustration of our model. Random switch- 
ing between states rji determines the current T)i<t) that interact 
with f(t) to change ice volume. The ice-sheet grows (b — 1) 
until it crosses a critical ice volume £, where it breaks up 
(b = —3). Once the ice-sheet is totally melted (V = 0.01£) 
growth starts again, (b) An example of a time series gener- 
ated by our model; note that the y-axis is inverted to allow 
easier comparison with Fig. [j]. The dashed lines indicate the 
maximal (f) and minimal (0.01£) ice-volume. 
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FIG. 4. (a) As for the Vostok record (Fig. |a), F(n) in- 
dicates random walk behavior (■). F(n) remains unchanged 
after the surrogate data test (gray triangles), (b) F(n) for 
the magnitude series reproduces the correlations observed for 
natural data (Fig. ^b). After the surrogate data test the mag- 
nitude series becomes uncorrelated indicating nonlinearity of 
our model, (c) The exponents r(q) versus the moment q for 
series generated by the model (10 realizations of 410kyr each, 
the average ± 1 standard deviation is shown) before (•) and 
after (o) the surrogate data test. The exponents measured for 
window scales between 0.8kyr and 25.6kyr. The nonlinear de- 
pendence of r(q) on q is similar to the natural data (Fig. ^c) 
and becomes linear after the surrogate data test, (d) The 
multifractal spectrum D(h) versus the exponent h(q). As in 
the data (Fig. ^|d) also here D(h) is broad before applying 
the surrogate data test and becomes narrower afterward. 
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